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We propose a multireference linearized coupled cluster theory using matrix product states (MPS- 
LCC) which provides remarkably accurate ground-state energies, at a computational cost that has 
the same scaling as multireference configuration interaction singles and doubles (MRCISD), for 
a wide variety of electronic Hamiltonians. These range from first-row dimers at equilibrium and 
stretched geometries, to highly multireference systems such as the chromium dimer and lattice 
models such as periodic two-dimensional 1-band and 3-band Hubbard models. The MPS-LCC 
theory shows a speed up of several orders of magnitude over the usual DMRG algorithm while 
delivering energies in excellent agreement with converged DMRG calculations. Also, in all the 
benchmark calculations presented here MPS-LCG outperformed the commonly used multi-reference 
quantum chemistry methods in some cases giving energies in excess of an order of magnitude more 
accurate. As a size-extensive method that can treat large active spaces, MPS-LCG opens up the 
use of multireference quantum chemical techniques in strongly-correlated ab-initio Hamiltonians, 
including two and three-dimensional solids. 

Introduction 

One of the most pressing theoretical questions in elec¬ 
tronic structure theory is how to deal with realistic 
strongly correlated electronic systems, which typically 
exhibit combinatorial complexity in the description of 
the ground-state wavefunction. There are two dominant 
paradigms in electronic structure algorithms, namely 
variational and projective methods, which for different 
reasons struggle to capture the entirety of the problem. 

Variational methods (such as and DMRGf^i^) are 
in principle robust techniques, but generally scale expo¬ 
nentially in the number of correlating orbitals and do 
not provide size-extensive energies (except in unreach¬ 
able exact limits), whilst projective methods, such as 
many-body perturbation theory and coupled cluster the¬ 
ory, spectacularly successful in describing weak corre¬ 
lation, fail as the underlying single-reference wavefunc¬ 
tion upon which they are build diminishes in impor¬ 
tance with growing strength of correlation. It is nat¬ 
ural to ask if a judicious combination of a variational 
and a projective method exists which can tractably han¬ 
dle realistic strong-correlation systems. Here we pro¬ 
pose such a method that is based on the linearized 
Coupled-Cluster (LCC) method^ and is implemented us¬ 
ing the matrix-product states (MPS)^ formalism to cap¬ 
ture highly multi-configurational zeroth- (and higher) or¬ 
der wavefunctions. 

To give a consolidated presentation, we first start by 
deriving the governing equations for the single reference 
LCC method which uses the Hartree-Fock wavefunction 
as the zeroth order state. We show that LCC is closely re¬ 
lated to the commonly used variational method, the con¬ 
figuration interaction with singles and doubles (CISD). 


We highlight the strengths and weaknesses of LCC com¬ 
pared to CISD. The main shortcoming of LCC (it be¬ 
comes divergent when near-degeneracies are present) can 
be overcome by using a multireference zeroth order wave- 
function. These multireference LCC equations have the 
same relation to the variational equations usually solved 
using DMRG^“— , as the LCC equations have to CISD. 
We show how the Matrix Product States can be used to 
efficiently solve the multireference equations by a small 
modification of the DMRG algorithm, resulting in a 
method which we call MPS-LCC. The resulting method 
is extremely powerful and we demonstrate its strength by 
solving several tough paradigmatic benchmark problems: 
first-row dimers at equilibrium and stretched geometries, 
the 1-band and 3-band (cuprate-like) Hubbard models in 
the strong-correlation regime C//t = 4 — 10 and the Cr 2 
dimer. 


Theory 

To derive linearized coupled cluster equations one 
starts with the coupled cluster singles and doubles wave- 
function written using the exponential ansatz |'I') = 
e^|4>o), where T = Ti -|- T 2 is the sum of the single and 
double excitation operators and |$o) is the Hatree-Fock 
wavefunction. When the CC wavefunction is substituted 
into the Schroedinger equation iJltk) = £'|4') and is left 
multiplied by e~^ we obtain e“^iLe^|$o) = A|4>o). Left 
projecting onto the Hartree Fock and a set of singly and 
doubly excited determinants |4>^) we obtain the expres¬ 
sion for the coupled cluster energy and the governing 
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equations for the t—amplitudes. 

E={%\e-^He^\^o) (1) 

0 =($Je-^i?e^|$o) (2) 

The set of non-linear Equations [2] can be solved to evalu¬ 
ate the t—amplitudes which can then be substituted into 
Equation [T] to obtain the coupled cluster energy. To ob¬ 
tain LCC equations the above equations can further be 
simplified by expanding the exponential using the Baker- 
Campbell-Hausdorff expansion and truncating the series 
at the first order, i.e e~^He^ = H + [H,T], to yield 

£;= ($o|i?|$o) + ($o|^|^'i) (3) 

0=($^|J?|$o) + ($mI(^-^^o)|^i) (4) 

Equations [3] and Equations |4] are the governing equa¬ 
tions of the single reference linearized coupled cluster the¬ 
ory, where we have defined |4'i) = T|$o) is the LCC 
correction to the Hartee-Fock wavefunction consisting of 
only single and double excitations. Equation 0] is now a 
linear equation in the unknown | 'I'l) which can be solved 
and substituted into Equation 0] to calculated the LCC 
energy. 

Now let us recall that the variational principle can be 
written as a set of equations — £'|'k) = 0, where 

1$^) are the basis states used to expand the wavefunction 
l^/). When the variational principle is used to optimize 
the CISD wavefunction we obtain 

E={<i>o\H\^o) + {^o\H\^i) (5) 

0 = + ( 6 ) 

where we have used intermediate normalization 
(($o|4') = 1) and as before |4'i) is the correction 
to the Hartree-Fock wavefunction consisting of only 
singly and doubly excited determinants. These equations 
are remarkably similar to Equations 0] and 0] with the 
small modification that the zeroth order energy Eq 
in Equation 0] is replaced with the variational energy 
E in Equation 0] This seemingly small change has a 
significant implication that the LCC energies unlike 
the CISD energies are fully size-extensive. The LCC 
energies are nearly equal to the full CCSD energies for 
weakly correlated problems. We demonstrate this by 
showing the correlation energy of diamond in Table U 
with an ab-initio Hamiltonian on a 2 x 2 x 2 fc—point 
sampling resulting in 64 electrons in 64 orbital problem 
(see [l^ for more details). 

The size-extensivity unfortunately comes at the cost 
of variationality and more problematic is the fact that 
the equations are prone to divergences^ in cases of near 
degeneracy when other determinants besides the Hartree- 
Fock have energies similar to Eq. 

This shortcoming can be overcome by formulating a 
multireference LCC method in which the Eq and 'kg in 
Equation 0] is replaced by the energy and wavefunction 


TABLE I: The correlation energy (Ej,/electron) of diamond 
calculated using various theories with an ab-initio Hamilto¬ 
nian. The DMRG calculation is still very far from convergence 
even though it was performed using an MPS with a large vir¬ 
tual bond dimension of M = 7500. The LCC energies are of 
comparable accuracy to the CCSD and the FCIQMC ener¬ 
gies and were calculated using a small MPS with M = 500, 
representing a speed up of over 3 orders of magnitude over 
the DMRC calculation. (See the main text for algorithmic 
details.) 


FCIQMC 

MP2 

CCSD 

LCC 

DMRC 

-0.6190 

-0.5145 

-0.6134 

-0.6235 

-0.5477 


obtained by fully correlating a set of orbitals around the 
Fermi surface. Such a multireference LCC method has 
been formally derived in the past by Bartlett et alJ^d^ 
as well as by Finbi^iii. Here, we use Fink’s formulation 
in which the single reference LCC is written as a pertur¬ 
bation theory using an ingenious use of a zeroth order 
Hamiltonian. This perturbation theory is then straight 
forwardly extended to multireferences cases thus result¬ 
ing in a multireference LCC theory. The advantage of 
Fink’s formulation is that it not only reduces to the LCC 
equations at the first order but systematic higher order 
corrections can also be generated. 

Following Fink’s formulation we first start by divid¬ 
ing all the orbitals into an active (or correlating) set, in 
which the orbital occupancies can be 0,1,2, a core set 
in which the occupancies are constrained to be 2, and a 
virtual set with zero occupancy. Then, we partition the 
full Hamiltonian H, whose ground-state wavefunction T' 
we seek, in terms of number-preserving operators within 
each subset: 

H = +y^{ij\kl)a\a\aiak = Hq + U (7) 

ij ijkl 

Hq = ^ + X! (8) 

ij', ijkl] 

Ariex—O Auex—O 

The constraint Auex = 0 implies that the operators do 
not transfer electrons between the three sets of orbitals, 
and U contains all remaining terms. 

To begin with we will assume that the ground-state 
eigenfunction of Hq can be found (later this as¬ 
sumption will be relaxed using the projector approxima¬ 
tion). This zeroth order wavefunction (which will in gen¬ 
eral have a combinatorial complexity) has an eigenvalue 
equal to the expectation value of the full Hamiltonian 
F;(o) = |iL and the first order energy is zero. It 
is possible to develop the usual perturbation theory mas¬ 
ter equations to express successive corrections (<&(")), to 
the wavefunction: T* = -1-4)^^)We expect 

this series to converge if the norms at each subsequent 
order rapidly diminish. In the limit in which the active 
space spans all orbitals, we recover full Cl (which is al¬ 
ways convergent), whereas in the opposite limit where 
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there are zero orbitals in the active space, we recover lin¬ 
earized Coupled Cluster theory, itself an excellent weak- 
correlation theory, but which diverges in strong correla¬ 
tion systems. We expect to have convergent theory for 
any level of correlation as long as the active space is suf¬ 
ficiently large. 

The equations governing are shown in EquationlHl 
where P is the projector on to the zeroth order wavefunc- 
tion {P = and Q is its complement (1 — P); 

the set of linear equations must be solved one at a time 
to obtain the order correction to the wavefunction 
Once is known, 2n and 2n -|- 1 order cor¬ 
rections to the energy calculated 

due to Wigners (2n -|- 1) rule using Equation fTOl 

(ffo = Q -f 

(9) 

n n—1 

£;(2ra) ^ _ EE 

fc=l j=l 
n n 

j^i2n+l) ^ _ EE 

fc=l j=l 

( 10 ) 

We use the variational principle for the perturbation 
theoryiSi^ which states that to solve Equation [9] it is 
sufficient to minimize the Hylleraas functional shown in 
Equation [T4l with respect to Only a small mod¬ 

ification to the DMRG algorithm, implemented in the 
Block code^ir— , was required to minimize the Hyller¬ 
aas functional. The details of the algorithm are outlined 
in the next sectional. 

Projector approximation: The cost of optimizing the 
zeroth order wavefunction scales exponentially with the 
size of the active space. This exponential cost can be 
circumvented by only approximately diagonalizing the 
zeroth order hamiltonian. The error in the zeroth order 
wavefunction can then also be corrected perturbatively. 
To do this we define a new zeroth-order hamiltonian 
Hq = PHqP + QHqQ, and the perturbing hamilto¬ 
nian becomes U = PHqQ -|- QHqP + U. It can then 
be shown that the use of the modified zeroth order 
hamiltonian changes the Equations [5] only at the first 
order to [Hq — where H is 

the unperturbed hamiltonian; and the expression of 
the second order energy changes from to 

($(0)!$(!)These results are remarkable because in 
essence they imply that even approximate solution of 
the zeroth order equation is sufficient and the original 
perturbation series can be used with only minor modifi¬ 
cations. 

Before moving on to describe our implementation, we 
would like to point out that Fink’s formulation is differ¬ 
ent than Bartlett’s formulation of LMRCC theory. Un¬ 
like Bartlett’s equations, in Fink’s equations the states 


with different number of electrons in the active, core and 
virtual spaces don’t interact with each other through 
the zeroth order Hamiltonian. Also, unlike Bartlett’s 
equations the present LMRCC equations cannot be ob¬ 
tained straightforwardly from governing equations of ic- 
MRC C^^i^^ by discarding terms that are higher order 
than linear in the excitation operator T. A key differ¬ 
ence between the current approach and the commonly 
used approach in multi reference methods is that we use 
a single MRS to represent the first order wavefunction 
instead of expanding it in the space formed by inter¬ 
nally contracted singles and doubles excitations out of 
the reference wavefunction. The states forming this space 
are not mutually orthogonal and without special care the 
equations can become ill-conditioned due to linear depen¬ 
dencies. The problem of linear-dependencies never arise 
in our formulation, further unlike the usual formulations 
at convergence our wavefunction is allowed to relax in the 
full uncontracted space of singles and doubles excitations. 
Another important difference is that unlike the usual for¬ 
mulation we do not need any reduced density matrices, 
for example up to sixth-order RDMs are required for a 
naive implementation of ic-MRCI calculations, although 
in practice this requirement can be realxed by using con¬ 
figuration state functions instead of internally contracted 
states22i^. The two other approaches used for incorpo¬ 
rating post-DMRG dynamical correlation are the Canon¬ 
ical Transformation theory and the methods based on ex¬ 
ploring the tangent space of the MPS^^— . The canon¬ 
ical transformation theory^li^^ tries to perform the uni¬ 
tary multi reference coupled cluster theory, but with the 
simplification that all RDMs higher than the two-body 
RDMs are evaluated using the cumulant approximation. 
The tangent space based methods are in principle quite 
similar to our current method with the main difference 
that the corrections to the reference MPS are restricted 
to linear combination of its tangent space vectors. This is 
a far more restrictive space than the one used here which 
is the one spanned by a single MPS with an arbitrarily 
large bond dimension. 

We would also like to emphasize that this method can 
be extended in several ways. First, the static one-particle 
and two-particle correlation functions of the ground 
state can be easily calculated^^i^. Second, in addition 
to the ground state, a set of low-lying excited states 
can also be calculated with a computational cost that 
scales linearly with the number of states using quasi¬ 
degenerate perturbation theory^^— . Finally, dynamical 
correlation functions can be calculated by combining the 
working equations of coupled cluster Green’s function 
framework^ with the dynamical DMR G^^d° . 


Implementation 

A general MPS representing a wavefunction jT) is 
shown in Eq. (HU, where Ui is the occupation of orbital i, 
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FIG. 1: A matrix product state (MPS) can be represented 
graphically using a series of 3-dimensional tensors in which, 
the physical index (pointing upwards) of the tensors denotes 
the occupation of the orbital and the other two indices, known 
as virtual indices, are sequentially contracted. Similarly, a 
matrix product operator (MPO) can be represented graphi¬ 
cally using a series of 4-dimensional tensors, with two physical 
indices and two virtual indices. The virtual indices of adjacent 
tensors are contracted sequentially. 



(Wl*) 





FIG. 2: The figure shown the overlap and a transition matrix 
element of an MPO Ho between MPS 4' and 4>. These are 
calculated by contracting the free physical dimension of the 
MPS and MPO sequentially as shown in the figure. By ap¬ 
propriately ordering the sequence of these contractions it can 
be shown that the cost of evaluating an overlap is 0{kM^) 
and a transition matrix element is O(k^M^) respectively. 


and ii ■ ■ ■ ik-i are the virtual indices that are contracted 
to obtain the final state. By increasing the size of the vir¬ 
tual indices an MPS can be used to represent any wave- 
function arbitrarily accurately. Similarly, operator O can 
be written in a matrix product operator (MPO) form as 
shown in Eq. (1131) . It should be noted that the expression 
in Equation [13] is very general, but it can be simplified if 
one limits the operator to have at most two body interac¬ 


tions. In such cases it can be shown that the virtual bond 
dimension of the operator need never be greater than k^, 
where k is the number of orbitals. The algebraic notation 
very quickly becomes cumbersome due to the rapid pro¬ 
liferation of indices; instead the graphical notation which 
is briefly explained below is much more convenient and 
intuitive. Here we can only give a short introduction to 
the notation, for more details we refer the reader to the 
excellent review article by Schollwocb^. 

Both the MPS and MPO can be conveniently repre¬ 
sented graphically as shown in Figure |T| In the figures 
each matrix (strictly speaking this is a three dimensional 
tensor) of the MPS is represented by a circle with three 
bonds jutting out, one pointing in the upward direction 
which corresponds to the physical index rii and two oth¬ 
ers pointing horizontally that correspond to the virtual 
indices. The bonds corresponding to the virtual indices 
of the adjacent matrices are joined together, which alge¬ 
braically corresponds to contracting the virtual indices, 
to obtain the wavefunction. The different MPS and MPO 
when written graphically are distinguished by the sym¬ 
bols used to represent their matrices; e.g. here we use a 
circle of Ik, a triangle for $ and a square for Hg respec¬ 
tively. 

It can be shown that taking the overlap between two 
MPS and calculating the matrix element of an MPO be¬ 
tween two MPS can be performed with a polynomial cpu 
cost of 0(kM^) and 0{k^M^) respectively (see Figure|l|). 
To get this computational scaling one needs to perform 
the various tensor contractions in a specific well-defined 
order. A suboptimal order of contractions can lead to 
computational cost that scales exponentially with the 
number of orbitals k. The partial derivative of overlap or 
operator expectation value with respect to one of the ma¬ 
trices of an MPS gives rise to a tensor which has exactly 
the same dimension as that of the matrix. This partial 
derivative in graphical language is represented by graph 
of the overlap or the expectation with the corresponding 
matrix removed from it as shown in Figure |3| 

In both MPS-LCC and DMRG the functional being 
optimized is quadratic in the wavefunction of interest. 
In the case of DMRG it is the energy functional, 

F;[T] =($|ff|$) - ($|F|$) (11) 

and in the case of LGG it is the Hylleraas functional 
shown in Equation 1141 
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($( 0 )|f/|$("- ^^)-E ^(fc)(|^( 0 )|^(n-fc)^ 


fc=l 


(12) 

(13) 


(14) 



d{^\^)/dAlY, 



d{m\^)/dAlY, 


FIG. 3: The figure shows the partial derivative of the overlap 
and transition matrix element with respect to the local ten¬ 
sor AT of the MPS. Each of these graphs represents 4M^ 
terms corresponding to taking the partial derivative with re¬ 
spect to each element of the tensor A"' • . 


In MPS-LCC the wavefunction of interest is written as 
an MPS and is then evaluated by minimizing the Hyller- 
aas functionals using the sweep algorithm. The key ele¬ 
ment of the sweep algorithm is that at each sweep itera¬ 
tion I only one tensor A()' is optimized keeping all the 
others constant. Figure [3] shows the partial derivative 
of (4'|$) and (fllliJld)) with respect to the unknown ten¬ 
sor A”' of wavefunction if. The governing equations 
that need to be solved at each sweep iteration are ob¬ 
tained by taking the partial derivatives of Equations [14] 
and equating them to zero. This converts a compli¬ 
cated multi-linear problem into a linear algebra problem 
(a linear equation) with the elements of tensor 
as unknowns. Standard iterative algorithms like the, 
Jacobi-Davidson and Conjugate-Gradient methods can 
be used to solve the linear algebra problems. By increas¬ 
ing the virtual bond dimension of the MPS expressing 
1$^"^), the Equation [T4| can be minimized arbitrarily ac¬ 
curately. The cpu cost per sweep iteration for calculating 
the two most expensive terms and 


($(”) |^|(I)("-i)) on the right hand size of Equation [T4| 
are 0{k‘^M^) and 0 (fc^M^M„_i) -|- 0{k'^MnM^_i) re¬ 
spectively, where Mn is the virtual bond dimension of 
the MPS representing the state The entire algo¬ 

rithm is implemented in the Block code which includes 
the ability to treat several different symmetries including 
SU(2) and non-Abelian point group. 


Benchmarks 

Earlier we showed that MPS-LCC is more efficient 
that the variational DMRG algorithm for weakly corre¬ 
lated systems like the diamond crystal. Here we demon¬ 
strate that it shows equally impressive performances for: 
first row dimers at equilibrium and stretched geometries, 
strongly correlated systems like 2 -dimensional I-band 
and 3-band Hubbard model at half filling and the Cr 2 
dimer. The half filled I-band Hubbard model is chosen 
because reliable APQMC results are available due to ab¬ 
sence of sign-problem. It should be pointed out that from 
the perspective of MPS-LCC calculations, half-filling rep¬ 
resents the hardest case and we expect the quality of re¬ 
sults to be better away from half filling. 


First row dimers 

We start by calculating the energy of the ground state 
of the C 2 dimer at various bond lengths using MPS-LCC 
with the double-zeta basis set. The MPS-LCC calcula¬ 
tions are performed on an MCSCE reference wavefunc¬ 
tion with an active space of ( 80 , 8 e) and the resulting 
energies are tabulated in Table [TTl The table also shows 
the PCI energies which are calculated by correlating all 
12 electrons in 28 orbitals using the DMRG algorithm as 
implemented in the Block code.The errors in the MC¬ 
SCE and MPS-LCC energies, calculated relative to the 
ECI energies, are plotted in EigurejH We see that there 
is a discontinuity in the MCSCE and MPS-LCC energies 
at a bond length of 3.10 bohr. This is because the 
and ^Ag energy curves cross between 3.05 bohr and 3.10 
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FIG. 4: The figure shows the error in the energies calculated 
using MCSCF and MPS-LCC methods relative to the FCI 
energies for the carbon dimer at various bond lengths using 
the cc-pVDZ basis set. There is a discontinuity in the energies 
of the MPS-LCC and MCSCF methods at 3.1 bohr because 
of a curve crossing between and ^Ag states. This curve 
crossing is shown in the inset. 


bohr bond lengths, with the former being the ground 
state at shorter bond lengths and the latter being the 
ground state at larger bond lengths. In our calculations 
the discontinuity in the MCSCF and MPS-LCC curves 
arise because we have only used the I? 2 /t subgroup of 
the full Dooh point group of the molecule. Besides the 
discontinuity at the curve crossing the MPS-LCC energy 
is both continuous and smooth despite the fact that at 
bond lengths greater than 3.05 bohr the ground state 
and the first excited states are nearly degenerate with a 
maximum separation of less than 6 mE^. 

We also benchmark the MPS-LCC method for C 2 , N 2 
and F 2 molecules at their equilibrium bond lengths of 
1.24253, 1.0977 and 1.4119 A respectively against the 
full configuration interaction (FCI) energies calculated 
using the FCIQMC method with up to quadruple-zeta 
basis set^. Here, the energies are calculated using vari¬ 
ous commonly used active-space methods like MRCI^i^, 
CASPT2, CASPT3ii, NEVPT3i^i^ and the MPS-LCC 
method developed in this work. For all these meth¬ 
ods two sets of calculations were performed, the first in 
which a complete active space configuration interaction 
(CAS-CI) wavefunction was used as a reference and in 
the second multi-configuration self consistent field (MC¬ 
SCF) wavefunction was used. Both CAS-CI and MCSCF 
calculations were performed with frozen core and all the 
valence orbitals included in the active space. The results 
of the calculations are shown in Table m and are plot¬ 
ted in Figure [5l These calculations show that MPS-LCC 
gives higher accuracy for these molecules compared to 
all other methods. One striking feature of these results 
is the fact that the quality of the MPS-LCC results are 


TABLE II: The table shows the ground state energy of the 
C 2 dimer calculated using various methods with the cc-pVDZ 
basis set. The MCSCF and MPS-LCC use an ( 8 e, 80 ) active 
space and the FCI energy is calculated by fully correlating 12 
electrons in 28 orbitals. 


r/ao 

FCI 

Energy/E^ 

MCSCF 

MPS-LCC 

1.80 

-75.4549 

-75.3501 

-75.4532 

1.85 

-75.5132 

-75.4080 

-75.5116 

1.90 

-75.5621 

-75.4564 

-75.5605 

1.95 

-75.6026 

-75.4965 

-75.6010 

2.00 

-75.6358 

-75.5294 

-75.6343 

2.05 

-75.6628 

-75.5560 

-75.6613 

2.10 

-75.6843 

-75.5771 

-75.6828 

2.15 

-75.7010 

-75.5935 

-75.6996 

2.20 

-75.7136 

-75.6058 

-75.7122 

2.25 

-75.7227 

-75.6145 

-75.7213 

2.30 

-75.7287 

-75.6202 

-75.7272 

2.35 

-75.7320 

-75.6232 

-75.7306 

2.40 

-75.7332 

-75.6240 

-75.7317 

2.45 

-75.7324 

-75.6229 

-75.7309 

2.50 

-75.7300 

-75.6202 

-75.7285 

2.55 

-75.7263 

-75.6161 

-75.7247 

2.60 

-75.7215 

-75.6109 

-75.7199 

2.65 

-75.7159 

-75.6047 

-75.7141 

2.70 

-75.7095 

-75.5979 

-75.7076 

2.75 

-75.7026 

-75.5904 

-75.7006 

2.80 

-75.6953 

-75.5825 

-75.6931 

2.85 

-75.6878 

-75.5743 

-75.6853 

2.90 

-75.6802 

-75.5659 

-75.6773 

2.95 

-75.6726 

-75.5573 

-75.6693 

3.00 

-75.6652 

-75.5487 

-75.6612 

3.05 

-75.6580 

-75.5400 

-75.6532 

3.10 

-75.6512 

-75.5208 

-75.6509 

3.15 

-75.6467 

-75.5163 

-75.6464 

3.20 

-75.6420 

-75.5116 

-75.6418 

3.25 

-75.6373 

-75.5068 

-75.6371 

3.30 

-75.6326 

-75.5020 

-75.6323 

3.35 

-75.6278 

-75.4972 

-75.6276 

3.40 

-75.6230 

-75.4924 

-75.6228 

3.45 

-75.6183 

-75.4877 

-75.6182 

3.50 

-75.6137 

-75.4831 

-75.6135 

3.55 

-75.6091 

-75.4785 

-75.6090 

3.60 

-75.6047 

-75.4740 

-75.6046 

3.65 

-75.6003 

-75.4696 

-75.6002 

3.70 

-75.5961 

-75.4654 

-75.5960 

3.75 

-75.5920 

-75.4613 

-75.5920 

3.80 

-75.5880 

-75.4573 

-75.5880 

3.85 

-75.5842 

-75.4535 

-75.5842 

3.90 

-75.5805 

-75.4498 

-75.5805 

3.95 

-75.5769 

-75.4463 

-75.5770 

4.00 

-75.5735 

-75.4429 

-75.5737 

4.05 

-75.5703 

-75.4398 

-75.5705 

4.10 

-75.5672 

-75.4367 

-75.5674 

4.15 

-75.5642 

-75.4339 

-75.5645 
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almost unaffected by the reference wavefunction. This 
is a well known feature of the CCSD method, but from 
these results it looks like the linearized version of the the¬ 
ory also shows this feature as long as the active space is 
large enough to avoid divergences. 

It can be seen that the results of the MRCI calculation 
for the F 2 dimer are much less accurate than the C 2 
and N 2 dimers. This is most likely due to the relatively 
small size of the zeroth order wavefunction which only 
has about 64 determinants in the active space compared 
to 4900 and 3136 determinants respectively for C 2 and 
N 2 respectively. The perturbation theories are in general 
somewhat less sensitive to the size of the Hilbert space 
of the zeroth order wavefunction because they are not 
variational. 

We have also performed calculations on the C 2 dimer 
with cc-pVQZ basis set at various bond lengths. Here 
the accurate benchmark data is obtained using the frozen 
core DMRG calculations^^ published recently. These re¬ 
sults again show the higher accuracy obtained by the 
MPS-LCC method relative to other methods. There are 
larger non-parallality errors at bond length of 1.6 A pos¬ 
sibly because of curve crosing between two Ag states (the 
intersecting and states belong to the Ag irre¬ 
ducible representation in the D 2 h subgroup) near this 
geometry. Quasi-degenerate perturbation theory can in 
principle be used to ameliorate these problems. 


Cr2 dimer 


The chromium dimer has been a challenging problem 
for quantum chemistry and large active spaces and basis 
sets are required to obtain the correct binding curve^^— . 
Here, we do not try to calculate the best binding curve 
that we can, but instead use some smaller benchmark cal¬ 
culations to compare the commonly used quantum chem¬ 
ical methods against the MPS-LCC method. In particu¬ 
lar, we carry out an all electron (48e, 42o) calculation on 
the Cr 2 dimer with an SVP basis set at a bond length of 
1.5 A. 

Table m shows the total energy and the well depth 
calculated using various quantum chemical methods, in¬ 
cluding coupled cluster with up to fourth order excitation 
and the contracted MRCI method^l as implemented in 
Molpro^. All the multireference calculations including 
the LCC were performed with a zero order wavefunction 
obtained by performing a (12e,12o) CASSCF calculation. 
The first order MPS-LCC wavefunction was represented 
with an MPS of virtual bond dimension 4000. It can 
be seen that the MPS-LCC method not only gives a to¬ 
tal energy closest to our best guess of FCI energy but 
the error in the calculated well depth is over an order of 
magnitude more accurate than any of the other methods 
shown here. 


u 


a 

u 


w 


o 

M 



FIG. 5: The error in the energies calculated using various ac¬ 


tive space methods relative to the highly accurate FCIQMC 
energies^ of C 2 , N 2 and F 2 molecules at bond lengths of 
1.24253, 1.0977 and 1.4119 A respectively. The upper panel 
shows the errors when the CAS-CI wavefunction was the ref¬ 


erence and the bottom panel used an MCSCF wavefunction 
as reference. In both cases the active space chosen was the 
full valence space containing eight orbitals. 
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FIG. 6: The error in the energies of the Carbon dimer using 
various active space methods relative to the energies calcu¬ 
lated using DMRG which themselves have an error of less 
than 0.1 mH relative to the FCI energies^. 






















TABLE III: The third column show the ground state energy and estimated uncertainty in Hartrees {Eh) of the C2, N2 and 
F2 molecules at bond lengths of 1.24253, 1.0977 and 1.4119 A respectively calculated using the FCIQMO^ method. The rest 
of the table shows the errors of various active space methods relative to FCIQMC in milli-Hartrees (mEh)- The active space 
used for these calculations consisted of the eight valence orbitals including 2s and 2p orbitals. Two sets of calculations were 
performed, one with the CAS-CI reference and the other with MCSCF reference. 


Molecule Basis FCIQMC 

CAS-CI reference (mEh) 

MCSCF reference (mEh) 

(Eh) 

MPS-LCC CASPT2 CASPT3 MRCI 

MPS-LCC CASPT2 CASPT3 MRCI NEVPT2 


C2 

dz 

-75.7285(1) 

0.0 

64.7 

30.8 

63.2 

1.4 

9.0 

5.2 

3.7 

21.4 

C2 

tz 

-75.7850(1) 

1.5 

73.7 

28.4 

65.7 

2.4 

10.6 

7.9 

7.2 

27.4 

C2 

qz 

-75.8023(3) 

1.9 

75.7 

25.1 

65.9 

2.6 

9.1 

8.5 

8.0 

27.2 

N2 

dz 

-109.2767(1) 

1.3 

34.3 

16.3 

20.7 

1.4 

18.4 

5.0 

7.0 

28.7 

N2 

tz 

-109.3754(1) 

3.7 

45.9 

18.2 

26.0 

2.6 

22.1 

7.7 

13.9 

37.1 

N2 

qz 

-109.4058(1) 

4.5 

49.0 

16.0 

26.6 

2.9 

20.7 

8.2 

15.9 

37.2 

F2 

dz 

-199.0994(1) 

3.1 

37.1 

19.8 

17.8 

4.3 

14.2 

8.4 

19.4 

44.2 

F2 

tz 

-199.2977(1) 

6.3 

56.6 

23.2 

25.0 

7.7 

19.2 

14.8 

35.4 

57.5 

F2 

qz 

-199.3598(2) 

6.7 

62.7 

22.0 

26.0 

7.9 

18.6 

16.0 

40.6 

59.6 


TABLE V: The table presents the absolute energies and the 
well depths (in E;,) calculated for the chromium dimer at 
1.5 A bond length with an SVP basis set. CASSCF, MRCI 
and MPS-LCC theories used an active space of 12 electrons 
in 12 orbitals. Notice that the well-depth calculated using 
the MPS-LCC method is over an order of magnitude more 
accurate than the result of any other method presented. The 
“FCI” energy was calculated by extrapolating a large DMRG 
calculation^ to zero discarded weight limit and is estimated 


r 

A 

DMRG 

(Eh) 

MPS-LCC MRCI CASPT3 CASPT2 NEVPTg 
(raEh) 

have a residual error 

of about 2 mEh- 


1.1 

-75.7613 

2.7 

8.2 

7.7 

9.7 

27.9 

Method 

1.5A(SVP) 


1.2 

-75.7992 

2.8 

8.3 

8.5 

9.5 

27.7 


Energy 

well depth 

1.24253 

-75.8027 

2.3 

8.4 

8.9 

9.5 

27.6 

CASSCF 

-2,086.2256 

0.167 

1.3 

-75.7994 

3.2 

8.5 

9.5 

9.5 

27.6 

CCSD 

-2,086.3880 

0.177 

1.4 

-75.7797 

3.9 

9.1 

10.9 

bo 

28.0 

CCSD(T) 

-2,086.4222 

0.150 

1.6 

-75.7241 

0.1 

12.4 

5.7 

22.8 

27.2 

CCSDTQ 

-2,086.4302 

0.143 

2 

-75.6460 

1.6 

8.4 

7.6 

15.4 

29.0 

MRCIC 

-2,086.4280 

0.138 


MPS-LCC 

-2,086.4349 

0.129 

FCI 

-2,086.4448 ± 0.002 

0.129 


TABLE IV: The second row of the table shows the DMRG 
energy in Hartrees for Garbon dimer at various bond lengths. 
The rest of the table shows the error in milli-Hartree (mi5/i)of 
various methods relative to the DMRG energies. All calcu¬ 
lations besides DMRG were performed on a reference wave- 
function obtained by a frozen core MGSGF calculation with 
eight electrons in eight orbitals active space. 


2d Hubbard model 

We first calculate the ground state energy of a half 
filled 18-site 2D Hubbard lattice at U/t = 4.0. In this 
system when 8 orbitals (4 degenerate orbitals above and 
4 below the Fermi surface in fc—space) are treated vari- 
ationally and the rest of the orbitals are correlated us¬ 
ing the MPS-LCC framework we get excellent agreement 
with the AFQMC results which are expected to agree 
with the FCI results to all shown significant digits. The 
zeroth order wavefunction calculated variationally only 
captures about 14% of the correlation energy, but with 
the inclusion of the first order correction to the wavefunc¬ 
tion we account for 101% of the remaining correlation 
energy. Even though we don’t necessarily recommend 


performing higher order perturbation theory, for illus¬ 
tration purposes we show that subsequent higher order 
corrections calculated up to the 7*^ order show rapid con¬ 
verge towards the FCI energy as shown in Table IVII To 
access the cost of the method, we show in Figure [7] that 
the third order MPS-LCC energy rapidly converges to 
its final value with an MPS bond dimension (M) of only 
about 200. This is to be contrasted with the extremely 
slow convergence of DMRG algorithm with both local¬ 
ized and delocalized fc—space basis. Given that the cpu 
time of both the MPS-LCC and the DMRG algorithm 
scale as 0{M^) we see about three orders of magnitude 
improvement in the computational cost. 
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FIG. 7: The graph shows the energy error of three p-MPS- 
LCC calculations, (where the approximate zeroth order wave- 
function is represented with an MPS of bond dimension 100 
(blue), 20 (green) and 10 (red) respectively) relative to the 
FCI energy versus the MPS bond dimension of the first order 
wavefunction. We have also shown the energy error of the 
DMRG calculation using localized orbitals (black) and delo¬ 
calized k-space orbitals (cyan) respectively. The MPS-LGC 
method shows several orders of magnitude speed up over the 
corresponding DMRG calculations. 

Now we assess the performance of the projector ap¬ 
proximation^ where the zeroth order wavefunction is only 
approximately evaluated. Here again the perturbation 
theory shows rapid convergence towards the FCI energy. 
In particular, when the energy of is minimized by 
using an MPS with a virtual bond dimension of only 20 
the errors in energy compared to full MPS-LCC rapidly 
diminishes as shown in Table IVTl 

Two p—MPS-LCC calculations, one with an (8e, 8o) 
and other with (24e, 24o) active space were performed 
on a 2-D Hubbard model with 50 sites. The zeroth order 
wavefunction in the two cases only accounted for about 
1% and 14% of correlation energy respectively. But re¬ 
markably, the third order correction to the energy was 
able to capture 95% and 99% of the remaining correlation 
energy. The first order correction to the wavefunction in 
the two cases above were represented by an MPS of bond 
dimension 5000 and 20000 respectively. Based on the re¬ 
sults of the smaller 18 site Hubbard model we expect this 
perturbation theory to be rapidly convergent although it 
was not possible to carry out these calculations due to 
the high computational cost. 


3-band Hubbard model 

Recently Schwarz et ali^ have published FCIQMG^i^ 
results on an undoped 3-band {p — d) Hubbard model 
with 10 unit cells. Each unit cell containing Cu02 
is represented by three orbitals, one 3dx2-y2 centered 
on the Cu atom and a 2px and 2py orbital on the 
the Oxygen atoms displaced in the x—direction and 
p—directions relative to the Cu respectively. For the 
details of the Hamiltonian we refer the reader to the 
original publication, but we would like to note that the 
model has inter-site potential and is extremely strongly 
correlated with the on-site repulsion divided by nearest 
neighbor hopping U/t 8. We performed an MPS-LCC 
calculation, where 10 holes were exactly correlated in the 
10 lowest energy orbitals and the effect of the remaining 
20 orbitals was taken into account perturbatively with 
a first order wavefunction represented by an MPS of 
bond dimension M = 1500. Table IVHI shows that the 
resulting MPS-LCC energy has an astonishingly small 
error of less than 0.0002 eV/hole compared to a very 
large DMRG calculation with an M = 6000. 


Conclusion and Outlook 

In this paper we show that the governing equations of 
multireference LCC theory can be efficiently solved using 
MPS by slightly modifying the DMRG algorithm. The 
theory has been used to obtain highly accurate energies, 
with a fraction of the cost of a variational DMRG calcula¬ 
tion, of several benchmark problems: ab-initio Hamilto¬ 
nian of diamond, to Cr 2 dimer, two-dimensional 1-band 
and 3-band Hubbard models at half filling in the strongly 
correlated regime ofC//t = 4 — 8. More broadly, our 
work demonstrates new possibilities for efficiently access¬ 
ing the ground state wave functions of highly correlated 
materials like transition metal oxides fully ab-initio with¬ 
out recourse to approximate models. 
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TABLE VI: The table shows the results of a MPS-LCC and p—MPS-LCC calculations (Energy/electron in units of t) to various 
orders of perturbation theory. For the 50 site Hubbard model we have performed two MPS-LCC calculations, one with an 
active space of 8 electrons in 8 orbitals (p—MPS-LCC(8)) and the other with 24 electrons in 24 orbitals (p—MPS-LCC(24)). 
We have also shown the coupled cluster results up to the CCSDTQP level, calculated using the program MRC C^°'^^ for the 
18 site Hubbard model (the CC calculations didn’t converge for the 50 site Hubbard model). 

Order 18 Site 2D Hubbard 50 Site 2D Hubbard 

of Theory MPS-LCC p-MPS-LCC CC p-MPS-LCC(8) p-MPS-LCC(24) 


0 

-0.804 

-0.802 

-0.778 

-0.679 

-0.705 

2 

-0.949 

-0.948 

-0.959 

-0.873 

-0.867 

3 

-0.960 

-0.961 

-0.965 

-0.867 

-0.878 
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-0.959 

-0.960 

-0.958 

- 

- 

5 
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-0.958 

- 

- 

6 

-0.958 

-0.958 

- 

- 

- 

7 

-0.958 

-0.958 

- 

- 

- 

FCI 


-0.958“ 



-0.880*’ 


“ Obtained by exact diagonalization, Auxiliary-field Monte CarloS, which has no sign problem at half filling. 


TABLE VH: Table shows the calculated ground state energy 
for the 3-band Hubbard model. Here Eq is the zeroth order 
energy obtained by fully correlating 10 holes in the 10 lowest 
energy orbitals. It is remarkable that the relatively much 
cheaper MPS-LCC theory is accurate to 4 significant places 
compared to the expensive DMRG calculation. 
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MPS-LCC DMRG 
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1.3399 

1.5821 1.5819 
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